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We presented an efficient algorithm, fast adaptive flat-histogram ensemble (FAFE) , to estimate the 
density of states (DOS) and to enhance sampling in large systems. FAFE calculates the means of an 
arbitrary extensive variable U in generalized ensembles to form points on the curve /3 S (U) = aS J^ , 
the derivative of the logarithmic DOS. Unlike the popular Wang-Landau-like (WLL) methods, FAFE 
satisfies the detailed-balance condition through out the simulation and automatically generates non- 
uniform (Pi,Ui) data points to follow the real change rate of j3 B (U) in different U regions and in 
different systems. Combined with a U— compression transformation, FAFE reduces the required 
simulation steps from 0(N 3 ^ 2 ) in WLL to 0(N 1 ^ 2 ), where N is the system size. We demonstrate 
the efficiency of FAFE in Lennard- Jones liquids with several TV values. More importantly, we show 
its abilities in finding and identifying different macroscopic states including meta-stable states in 
phase co-existing regions. 
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Complex systems, such as proteins, glasses, and many 
materials in the phase-coexisting regions, can be charac- 
terized as having a rugged potential energy landscape. 
The molecular simulations of such systems face one cen- 
tral problem: how to generate samples that cover enough 
of the energy landscape to characterize the properties of 
the systems correctly. Normal Metropolis Monte Carlo 
and molecular dynamics techniques often only cover part 
of the conformational space because of the existence of 
high free-energy barriers. Over the last decade, some 
enhanced sampling techniques, such as the replica ex- 
change methods (REM) Jl|_and many flat-histogram 
techniques have been developed 

in attempts to systematically cover more conformation 
space. In particular, the flat-histogram techniques pre- 
scribe a collective variable and apply modified distribu- 
tion functions to flatten the visited histogram in this one- 
variable space. 

Here we propose the density of states (DOS) to gen- 
erate the flat-histogram ensemble. The DOS is defined 
as ft(U) = J dr5(U — U(r)), where r is the conforma- 
tion vector and U (r) can be any collective variable. We 
define S(U) = mft(J7), and p,(U) = While U 

is the potential energy, S(U) would correspond to the 
entropy and (3 S (U) is the inverted statistical tempera- 
ture. While the visited histogram of simulations in a 
S(U(r)) ensemble, which is defined as P(r) cx e~ s ( u ( r ^ 
H(U) oc e S(U)-s(if)^ ig amiost fl at) s(U) can be well es- 
timated by S(U) = S(U) + h\H{U). We can generate a 
sample in the S(U) ensemble, and calculate the mean of 
an arbitrary variable Q(r) in an arbitrary ensemble f(X) 
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k th conformation of the sample [13 |. 
weight factor defined as: = 
e S(u(r k ))-f(x(,r k ))_ The gtatistical error of tnis we i g hted 

mean is proportional to B -1 / 2 , where B is the effective 
size of the weighted sample, B = (J2 w k) 2 1 J2 w k- Only 
when sufficient conformations of the sample are located 
inside the important regions of the f(X) ensemble (i.e., 
B is large enough), will the weighted sample provide a 
good estimate of the ensemble mean. 

The Wang-Laudau (WL) algorithm Q estimates S(U) 
by updating S(U) from an initial guess until the vis- 
ited histogram is roughly a constant. The WL algo- 
rithm can be easily implemented and has been widely 
applied 0, 11, [l2|. However, two features of this algo- 
rithm would limit its applicability. First, the required 
number of simulation steps is at least 0(N 2 ) jjjj], where 
N is the size of system; hence simulations of large sys- 
tems will be costly. This estimate arise from the fact 
that the area below the curve S(U) is order- N 2 (because 
both S and U are proportional to N), which must be 
filled by using the small kernel functions with presumed 
width and height in WL Q. Second, detailed balance is 
not satisfied in the WL algorithm [6(, which leads to po- 
tential numerical instabilities because of the correlations 
between successive conformations in the simulations. In 
a recent extension to the original WL, the statistical- 
temperature molecular dynamics (STMD) method [§], 
the intensive variable T(U) = /3j 1 ([/), instead of the ex- 
tensive variable S(U), is updated at each simulation step. 
This method is faster than the original WL, since we do 
not need to fill the iV-order height of the curve S(U) due 
to the replacement. However, STMD still needs to di- 
vide the extensive variable U into order- N homogeneous 
bins with a prescribed bin size and it also does not sat- 
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isfy detailed balance. Thus the needed simulation steps 
are estimated to be 0(N 3 ^ 2 ). Another improvement of 
WL, the multimicrocanonical method (MMCM) [fjj], cal- 
culates the microcanonical temperature T m (E) from the 
microcanonical mean, then relates the reversed T m (E) to 
(3^ 1 (E) based on the thermodynamical equivalence of dif- 
ferent temperature definitions. The MMCM has shown 
better numerical stability, but is also order N 3 / 2 because 
it similarly needs an TV-order uniform grid in U before 
starting a simulation. Extending MMCM to the cases 
where the thermodynamical equivalence can not be used 
(e.g., U is not the total potential energy of system), is 
nontrivial. 

In this letter, we present a more efficient algorithm, 
fast adaptive flat-histogram ensemble (FAFE), to accu- 
rately estimate the DOS and enhance sampling in large 
general systems. FAFE calculates the ensemble means 
and fluctuations of U using a series of generalized en- 
sembles, then adaptively forms the points [/,;) on the 
curve P S (U). These Ui points non-uniformly distributed 
in U space and the intervals are dependent on both U and 
N to follow the real shape of the curve f3 s (U): we auto- 
matically form many U points in regions where f3 s (U) 
varies quickly, and few Ui points in the regions where 
f3 s (U) varies slowly. A compression transformation of 
U is also applied to reduce the visits to the U regions 
where (3 S (U) has been well estimated. Therefore the re- 
quired simulations steps of FAFE to estimate S(U) are 
not more than N 1 / 2 . Due to the good numerical stabil- 
ity of FAFE, we can visit and identify coexistent phases 
(including meta-stable phases) by constructing the DOS 
for each phase in the system, where the Wang-Landau 
like (WLL) methods cannot work unless the exact order 
parameters are known. 

In a f(U) ensemble, the histogram, fly (17) oc 
e s(u)-f(u) ^ can jj e approximated as the sum of a few 
Gaussian functions, Hf(U) oc ^ Ci exp{ — ^(U — U 1 ) 2 }, 
in the limit of large N. If Hf(U) is a single-peak function, 
the mean of U in the f(U) ensemble, (U)f, is approxi- 
mately equal to the maximal point of Hf(U), we have, 

-djjHu)f ~ / +Qfj2\(u) f > ( 2 ) 

where, aj is the fluctuation of U in the f(U) ensemble. 
Thus, we can calculate (U) / and ct/ in different f(U) en- 
sembles to generate a series of points on the curve S (U). 
The complete curve of (3 a (U) can be interpolated from 
these points, and S(U) can be estimated by integrating 
Ps(U). This approximation, eq.([2|), is sufficient in U re- 
gions where simulations based on the constructed S(U) 
can visit sufficiently. 

Different f(U) ensembles, such as f(U) — (3U (canon- 
ical ensembles), and = (E — U)~Q(E — U) (mi- 



crocanonical ensembles), can be used to construct /3 S (U). 
It is also possible to use generalized ensembles, such as 

f(U)={3 U+^a(U-U ) 2 , (3) 

to generate points on (3 S (U). Because the derivative 
is positive, phase-coexistence regions are unstable in 
canonical ensembles. But in a generalized ensemble, e.g., 
a > t^j-, they may be stabilized. Thus, once we choose 
suitable parameters for (/3q, Uo,a), we can even generate 
the j3 s (U) data points in the phase-coexitence regions. 

The implementation of FAFE is rather straightfor- 
ward. Without losing any generality, we choose U as 
the total potential energy, and define the statistical tem- 
perature T S (U) = P~ X (U) to illustrate it: (i) Choose an 
initial T and set T S (U) = T , then S(U) = U/T ; (ii) 
Simulate a segment of trajectory in the current S(U) en- 
semble and generate an equilibrium sample; (iii) Choose 
a few f(U) functions (e.g. f(U) = U/Ti for a few differ- 
ent Ti ) and calculate the ensemble means U and fluc- 
tuations <Ji from the generated sample at (ii) based on 
the eq.©. We eh oose Ti from the current T S (U) func- 
tion: the corresponding Ui, which satisfies the equation 
Ti = T S (U), must be located in the current visited re- 
gions. In addition, \Ui — J7j-i| ~ Uj— i, where er^-i is the 
fluctuation of U at the temperature Tj_i. Moreover, we 
only need to apply the mean and fluctuation at Ti when 
the effective size of the weighted sample, B, is sufficiently 
large. Note that we can calculate the means and fluctua- 
tions at many Ti from the same sample without additional 
simulations. The CPU time required for calculation of 
the weighted mean is almost negligible compared to that 
for generating samples; thus we can generate more data 
points without much additional computational cost; (iv) 
From the discrete points (Ui,Ti) (and the derivative of 
T S (U), T( = T 2 /a 2 ), interpolate T S (U) and extrap- 
olate it linearly beyond the last U [l5(. Integrate S(U) 
from the new T S (U); (v) Repeat (ii), (iii) and (iv) until 
the T S (U) and S(U) curves are complete. In Step (iii), 
we can check for the possible multi-peak characteristics 
of the weighted sample by using the generalized ensem- 
bles described by eq.([3]). For example, we first choose 
a = (i.e. canonical ensembles) to calculate (Tj,C/j), 
then set [3$ — Tr , Uq — U and use a positive a to cal- 
culate the means in the new f(U) ensemble, comparing 
the results to determine if larger a should be used to de- 
press the possible multi-peak effects. In practice, unless 
simulations are unexpectedly found to be trapped, this 
checking step is not necessary. The points (fii, Ui) can be 
calculated by using many samples generated in different 
segments of simulations in (ii) rather than only using a 

single segment of simulations, U <— , where £P 

and fl J are the the weighted mean of U and the corre- 
sponding B factor from the j th sample, respectively, and 
the statistical error of U, which is measured by the total 
B = & factor, can be depressed. 
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FIG. 1: The T S (U), c v = (NT^U))' 1 and S(U) in systems 
with different size N are calculated in a large temperature 
interval. The data points are adaptively generated to form 
the flat-histogram ensemble with the required accuracy, so 
that more data points are formed in larger systems and in the 
lower temperature region. 



In many cases, the T S (U) may be well estimated lo- 
cally, (i.e., the B factor of the corresponding U is large 
enough), before the other regions have been explored. 
For example, in usual physical systems, T S (U) at high 
temperature can be obtained before the low tempera- 
ture region is visited. The visit frequency in the large- 
U regions, which is proportional to N even in the flat- 
histogram ensemble, can be further depressed to force the 
simulation to focus more on the low-f region. To do so, 
we transform U to a new variable Y by ^ = g c ( U ~^ J " ), 
where the compression function g c (z) slowly approaches 
zero while z — > oo, but unity while z < 0. Here U c and b 
are parameters. One possible choice is g c (z) — e f ^ + ^ c 
with the parameter c ~ 2 or 3. At the end of each seg- 
ment of simulations, we reset U c at the lower bound- 
ary of the explored region and choose b to make the 
expected visit probability above U c comparable to that 
below U c , then use S(Y) to generate the flat-Y sample. 
Here S(Y(U)) = S(U)-\ng c (¥^), and the correspond- 
ing P(U) ~ 1/(17 - U c ) while U > U c . The flat-Y sim- 
ulations can still reach the high-?/ region to be ergodic, 
but mainly visit the low-C/ region to estimate T S (U) in 
the unexplored region. This compression transformation 
significantly decreases the computational cost of our al- 
gorithm for very large systems. 

To illustrate FAFE, we apply it to truncated-and- 
shifted Lennard-Jones (LJ) fluids (r c = 2.5ay) at re- 
duced density p — 0.8. We choose U as the total poten- 
tial energy, and aim to calculate DOS between T) = 0.5 
and Th = 100, which corresponds to a very large U range 
in large systems. The ensemble mean of U is sufficiently 
exact if the error in U is smaller than 5 percent of the 
fluctuation of U (i.e., B > 400). We apply molecular 
dynamics (MD) simulations with a Langevin thermostat 
and save a conformation every 10 MD steps. T S (U) is 
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FIG. 2: The efficiency of our algorithm: the required MD 
steps grow as N 1 ^ 2 . 




FIG. 3: T S (U) at super-saturated liquid and solid phases in 
LJ system with N = 108, respectively. The transition can 
easily happen while U < U\ (from liquid to solid) and while 
U > U2 (from solid to liquid). 



calculated and updated every 10 5 steps, when the pa- 
rameters of the compressive transformation are updated. 
In our simulations, the distance between neighboring [7j 
points is far larger than that in WLL [§] in high-T regions 
and large- N systems, but automatically decreases in low 
T regions and small TV systems. Figure (JTJ shows that 
T S (U/N), c sv = (dT s /d(U/N))- 1 and s(u) = S(U/N)/N 
of systems with several different N values are quickly 
generated in the large temperature regions. The curves 
for different N all collapse into single curves as expected. 
Figure © shows the efficiency of FAFE, the number of 
MD steps t oc N" 1 with 7 w 0.5. 

As shown in Fig. in the lower temperature region, 
T S (U) has multiple branches that correspond to the DOS 
of some separated parts in conformational space. For a 
segment of simulated trajectory, we may calculate the 
corresponding (T, (U)t) points, and learn which branch 
of T S (U) the points are located on (i.e., which macro- 
scopic state the segment of trajectory resides in). The 
trajectories generated under the high-J7 branch of T S (U) 
(liquid phase) have the expected flat histograms in the 
region U > U\, but quickly transition into the \cw-U 
region (solid phase). Similarly, simulations under the 
low-C/ branch of T S (U) (solid phase) visit the U < U2 
region with equal probabilities, but transition into the 
liquid phase. The overlap of the two branches indicates 



that the potential energy of the liquid phase can be lower 
than that of the solid phase, and f is definitely not the 
proper order parameter of the liquid/solid phase transi- 
tion. It is very hard to visit (and more difficult to iden- 
tify) the low-energy liquid and high-energy solid confor- 
mations in current existing simulation techniques. Thus, 
FAFE, which is able to visit more parts of the configura- 
tion space and identify them by forming their own DOS 
without requiring the prior knowledge of precise order 
parameters, might provide new promise in understand- 
ing phase transition processes. 

It is possible to calculate the total 5(f) and T S (U) 
from these obtained T a {U) branches [l6j]- However, even 
in the exact 5(f) ensemble, it is still not guaranteed that 
simulations can freely transition between the states. In 
this case, the visited probability of conformation inside 
each iso-f is constant, P{r) oc e~ s ^ u \ then any small- 
fraction part in an iso-f subspace is still rarely visited 
although the total visited probability of the iso-t/ sub- 
space is not small. If a small-fraction part corresponds 
to the passage between two states, simulations in the 
5(f) ensemble still do not easily pass the passage. The 
5(f) ensemble only provides chances to search transi- 
tion passages in more iso-t/ subspaces. This is a general 
limitation in all the low-dimension flat-histogram ensem- 
ble methods. The algorithms with good stabilities (e.g., 
simulations satisfying detailed balance) are helpful in an- 
alyzing and (partly) overcoming the limitation. For ex- 
ample, FAFE generates the T S (U) and 5(f) (with ar- 
bitrary constant) in each visited region instead of the 
expected total T S (U) and 5(f) of the whole conforma- 
tional space. Thus, even when the equilibration among 
the regions is difficult, it is still possible to provide many 
useful results. From the branches of 5(f) or T S (U), it is 
possible to construct the 5(f) for whole the conforma- 
tional spac In addition, we can apply the general- 
ized ensembles described by eq.Q to find more (possibly 
metastable) states. For example, besides the stable liq- 
uid and solid phases in the LJ system, we detect several 
metastable solid states with different crystalline struc- 
tures as well as the unstable multi-phase coexistence, 
which are shown in Fig.Q. These complex solid phases 
arise from the fact that the system is fixed at p = 0.8 so 
the interparticle distance in the face-center cubic struc- 
ture does not match the equilibrium distance of LJ po- 
tential, r = 2 1 / 6 . 

The flat-histogram ensembles based on multiple col- 
lective variables are expected to be better for enhanced 
sampling than single variable 0,113. The required com- 
putational cost will be 0(N ld ) 1 where d is the num- 
ber of the collective variables. FAFE (7 = 1/2), which 
is much faster than WLL (7 = 1.5 ~ 3), is efficient 
in sampling large complex systems, such as glasses and 
proteins. It is useful to construct the canonical DOS 
Q(U;T) = / e- v ^/ T 8(U-U(r))dr instead of the normal 
DOS in FAFE in order to keep the connectivity of macro- 
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FIG. 4: The T S (U) function in a LJ system with N = 108 
and p — 0.8 at the low temperature region. Each branch of 
the function corresponds a macroscopic phase. The unstable 
liquid/solid coexistence region, where T' S (U) = (iVci,) -1 are 
negative, is shown by the blue filled squares. The green down- 
triangles corresponds to the complete face center cubic solid 
which is a metastable state in the current p. Some phase tran- 
sitions from metastable solids to the stable solid (the black 
circles) can be found at lower temperature. While T„(U) — > 0, 
more complex phase behaviors are expected. 

molecules and more focus on the desired temperature T, 
if the total potential energy is not chosen as f . It is 
also possible to combine FAFE with the replica exchange 
method (REM): we generate the branches of T S (U) and 
use them in different replicas and exchange conforma- 
tions among these replicas; thus we can use O(N ) repli- 
cas in large systems instead of the N 1 / 2 — order replicas 
required in conventional REM, which is the main limita- 
tion of applying REM. 
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